! DEREK---comments

module mObject

        use mPrecision
        use mConstants
!        use mParameters

        implicit none

        integer, parameter :: MAX_ROWS = 7, &      ! Max # of rows used in Richardson extrapolation.
                              MAX_EQUATIONS = 2, & ! Max # of eqs in system of 1st order diff. eqs.
                              DIFFEQ_MAX_PARAMS = 3
        
        type tDiffEq
                real(D) :: alpha(MAX_ROWS, MAX_ROWS) ! alpha(k, q) := convergence correction factors, NR Eq. (16.4.10).
                ! A(i) := estimate of computational work needed to get to row k of the extrapolation tableau
                ! (cf. i.e., section 3.1 of Numerical Recipes in Fortran):
                real(D) :: A(MAX_ROWS + 1)
                real(D) :: param(DIFFEQ_MAX_PARAMS)
                real(D) :: y(MAX_EQUATIONS)
                real(D) :: x
                real(D) :: eps
                real(D) :: h
                real(D) :: hMin
                integer :: nEq   ! Number of 1st order differential equations ( <= MAX_EQUATIONS)
                integer :: k_max ! Optimal row number limit from NR Eq. (16.4.13).
                integer :: k_opt ! Optimal row number from NR Eq. (16.4.9).
                integer :: nGood, nBad ! Number of good (converged at attempted stepsize) and bad (required stepsize adjustment) steps taken.
                logical :: first
        end type tDiffEq

        type tUniverse
                real(D) :: omega_m0
                real(D) :: omega_v0
                real(D) :: w_0
                real(D) :: w_a
                real(D) :: omega_b0
                real(D) :: h
                real(D) :: A_s
                real(D) :: sigma_8
                real(D) :: b_1
                real(D) :: b_2
                real(D) :: growthFactorNorm
        end type tUniverse

        integer, parameter :: TRANSFER_FUNCTION_LEN = 2, &
                              NORMALIZATION_LEN = 2

        type tLPS
                type(tUniverse) :: universe
                real(D) :: mass_nu
                real(D) :: spectralIndex
                real(D) :: run
                real(D) :: norm
                character(len = TRANSFER_FUNCTION_LEN) :: transferFunction
                character(len = NORMALIZATION_LEN) :: normalization
        end type tLPS

        type tFitParams
                type(tLPS) :: lps
                type(tUniverse) :: universe
                real(D) :: param
        end type tFitParams

        integer, parameter :: MAX_NUMORDER = 5
        integer, parameter :: MAX_DENORDER = 5

        type tRFit
                real(D) :: coeffs(MAX_NUMORDER + MAX_DENORDER + 1)
                real(D) :: invCoeffs(MAX_NUMORDER + MAX_DENORDER + 1)
                real(D) :: ll, ul
                integer :: numorder
                integer :: denorder
                logical :: inverse
                logical :: xLog
                logical :: yLog
                logical :: fitDone
                type(tFitParams) :: params
        end type tRFit

        integer, parameter :: MAX_GAUSS_COEFFS = 50

        type tCFit
                real(D) :: coeffs(MAX_GAUSS_COEFFS)
                real(D) :: ll, ul
                integer :: numCoeffs
                logical :: xLog
                logical :: yLog
                logical :: fitDone
                type(tFitParams) :: params
        end type tCFit

        type tNu
                type(tRFit) :: nuFit
                type(tCFit) :: dNuFit
                type(tLPS) :: lps
!               Set this to .true. if nu0 depends on the scale factor, e.g. for
!               the Eisenstein & Hu transfer function and massive neutrinos:
                logical :: nu0ADependence
                real(D) :: a, dNuFit_a
        end type tNu

        type tHalo
                type(tNu) :: nuFit
                real(D) :: fracSub
                real(D) :: mu
                real(D) :: subLL, subUL
                real(D) :: biasNorm, normalizedForA
                real(D) :: m_star
                logical :: cDistribution
        end type tHalo

        type tGSV
                type(tLPS) :: lps
                type(tCFit) :: GSVFit
        end type tGSV

        type tSmith
                type(tGSV) :: GSV
                type(tCFit) :: dGSV
                type(tCFit) :: ddGSV
                real(D) :: a
                real(D) :: rmid, rneff, rncur, rknl, om_m, om_v
                logical :: useLinear
        end type tSmith

        type tChi
                type(tRFit) :: chiFit
                type(tUniverse) :: u
        end type tChi

        type tGalaxyBin
                type(tChi) :: chi
                type(tLPS) :: lps
                real(D) :: galaxyDensity
                real(D) :: fracSky
                real(D) :: z_0
                real(D) :: zMean
                real(D) :: zSigma
                real(D) :: zLL
                real(D) :: zUL
        end type tGalaxyBin

        integer, parameter :: INTEGRAL_CHAR_LEN = 2
        integer, parameter :: INTEGRAL_MAX_PARAMS = 5
        integer, parameter :: INTEGRAL_MAX_INTEGER_PARAMS = 3
        integer, parameter :: INTEGRAL_MAX_DIMENSION = 12

        type tIntegral
                type(tSmith) :: smith
                type(tHalo) :: halo
                type(tGalaxyBin) :: galaxyBin
                type(tChi) :: chi
                type(tLPS) :: lps
                type(tUniverse) :: u
                real(D) :: params(INTEGRAL_MAX_PARAMS)
                real(D) :: ll, ul
                real(D) :: eps
                integer :: dimension ! matrix integrand is dimension x dimension (internal storage)
                integer :: integerParams(INTEGRAL_MAX_INTEGER_PARAMS)
                logical :: converged_2(INTEGRAL_MAX_DIMENSION, INTEGRAL_MAX_DIMENSION)
                logical :: open
                logical :: converged
                logical :: error
                character(len = 3) :: changeVariables
                character(len = INTEGRAL_CHAR_LEN) :: char
        end type tIntegral

        ! Default initial spacing for derivative:
        integer, parameter :: DERIVATIVE_MAX_PARAMS = 4

        type tDerivative
                type(tGalaxyBin) :: galaxyBin
                type(tChi) :: chi
                type(tLPS) :: lps
                type(tUniverse) :: u
                real(D) :: params(DERIVATIVE_MAX_PARAMS)
                real(D) :: h
                real(D) :: error
        end type tDerivative

        integer, parameter :: ROOT_N_PARAMS = 2

        type tRoot
                type(tSmith) :: smith
                type(tChi) :: chi
                type(tLPS) :: lps
                real(D) :: params(ROOT_N_PARAMS)
                real(D) :: ll, ul
                real(D) :: tolerance
                real(D) :: output
                logical :: converged = .false.
        end type tRoot

        private
        public :: MAX_ROWS, MAX_EQUATIONS, DIFFEQ_MAX_PARAMS, tDiffEq
        public :: tUniverse
        public :: TRANSFER_FUNCTION_LEN, NORMALIZATION_LEN, tLPS
        public :: tFitParams
        public :: tRFit
        public :: MAX_GAUSS_COEFFS, tCFit
        public :: tNu
        public :: tHalo
        public :: INTEGRAL_MAX_PARAMS, INTEGRAL_MAX_INTEGER_PARAMS, INTEGRAL_CHAR_LEN
        public :: tIntegral
        public :: tChi
        public :: tGalaxyBin
        public :: tGSV
        public :: tSmith
        public :: tDerivative, DERIVATIVE_MAX_PARAMS
        public :: tRoot, ROOT_N_PARAMS

end module mObject

